function DecodingProbability_v27_dejan()

% Clear commandline
clc
close all

% Time it
tic

% Transmitter range
T_min=1;
T_step=1;
T_max=70;

% Field cardinality is the q vector input
% Layer selection probability is the g vector input

% Layer dimensions (these are dejan examples)
k1=15;
k2=25;
K1=k1;
K2=k1+k2;

% Gamma values that shall be plotted.
g1=[0.3 0.7];
g2=[0.4 0.6];
g3=[0.5 0.5];
g4=[0.6 0.4];
g5=[0.7 0.3];

q=2^8;

% Create all transition matrices
global PP_q1 PP_q2
PP_q1=cell(1,K2);
PP_q2=cell(1,K2);
for i=K1:K2
    PP_q1{i}=TransitionMatrix(i,q);
    PP_q2{i}=TransitionMatrix(i,q);
end

disp('Done with all transition matrices')

% Calculate decoding probabilities #set 1
l1_prob_1=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g1);
l2_prob_1=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g1);

% Calculate decoding probabilities #set 2
l1_prob_2=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g2);
l2_prob_2=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g2);

% Calculate decoding probabilities #set 3
l1_prob_3=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g3);
l2_prob_3=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g3);

% Calculate decoding probabilities #set 4
l1_prob_4=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g4);
l2_prob_4=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g4);

% Calculate decoding probabilities #set 5
l1_prob_5=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g5);
l2_prob_5=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g5);

% Plot nice graph!
hold all

% Linespec options here
% http://www.mathworks.se/help/techdoc/ref/linespec.html

plotter(l1_prob_1,l2_prob_1,T_min,T_step,T_max,K1,K2,g1);
plotter(l1_prob_2,l2_prob_2,T_min,T_step,T_max,K1,K2,g2);
plotter(l1_prob_3,l2_prob_3,T_min,T_step,T_max,K1,K2,g3);
plotter(l1_prob_4,l2_prob_4,T_min,T_step,T_max,K1,K2,g4);
legend_h=plotter(l1_prob_5,l2_prob_5,T_min,T_step,T_max,K1,K2,g5);

% Set proper legend
label_1=strcat('L1, \Gamma_1=',num2str(g1(1)));
label_2=strcat('L2, \Gamma_1=',num2str(g1(1)));
label_3=strcat('L1, \Gamma_1=',num2str(g2(1)));
label_4=strcat('L2, \Gamma_1=',num2str(g2(1)));
label_5=strcat('L1, \Gamma_1=',num2str(g3(1)));
label_6=strcat('L2, \Gamma_1=',num2str(g3(1)));
label_7=strcat('L1, \Gamma_1=',num2str(g4(1)));
label_8=strcat('L2, \Gamma_1=',num2str(g4(1)));
label_9=strcat('L1, \Gamma_1=',num2str(g5(1)));
label_10=strcat('L2, \Gamma=_1',num2str(g5(1)));

set(legend_h,'String',[label_1;label_2;label_3;label_4;label_5;label_6;label_7;label_8;label_9;label_10],'location','NorthWest')
set(legend_h,'interpreter','tex')
set(legend_h,'fontsize',8)

% Save plot
figname=strcat('dejan.eps');
print(gcf,'-deps',figname)

% Save data
save('uep_ew_analytic_dejan')

% Time it
toc

end

% Calculate layer 1 probability
function l1_prob = layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g)

l1_prob=zeros(T_max,1); % Decoding 1. Layer probabilities

for tx =T_min:T_step:T_max % For each number of recv packets
    
    sol=0;
    
    for n = 0:tx % For each permutation of a number of recv packets
        
        % Layer 1 by itself
        val=PM2(n,K1,K1,q);
        
        val_test=0;
        
        % Layer 1 and Layer 2 gives rank K2 (This way we also get L1!)
        % Sum prob for all possible ways to achieve rank K2 with given permutation of recv packets
        
        for i=0:K1-1
             
            tmp1=PM2(n,K1,i,q);   
            
            %Optimization We only need to calculate tmp2 if tmp1!=0
            if tmp1==0
                continue
            end
            
            tmp2=PM2(tx-n,K2-i,K2-i,q);
            val_test=val_test+tmp1*tmp2;
            
        end
        
        % We have counted all the ways L1 can become full rank by itself
        % We have counted all the ways L2 can become full rank (except when
        % L1 is full rank!)
        % Both outcomes are valid for getting L1 and since they are disjoint we can just add them!
        % P(a)+P(b)-P(ab), where P(ab)=0 because they are disjoint!
        
        sol=sol+(val+val_test)*binopdf(n,tx,g(1));
        
    end
    
    l1_prob(tx)=sol;
    disp(['Layer 1: ' num2str(tx) ' out of ' num2str(T_max)])
    
end

end

% Calculate layer 2 probability
function l2_prob = layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g)

% Decoding 2. Layer probabilities
l2_prob=zeros(T_max,1);

for tx =T_min:T_step:T_max % For each number of recv packets
    
    sol=0;
    
    for n = 0:tx % For all permutations of recv packets
        
        val=0;
        
        for i=0:K1 % For all possible ways to achieve rank K2 with given permutation of recv packets
            
            tmp1=PM2(n,K1,i,q);
            
            % Optimization no need for tmp2 when tmp1=0!
            if tmp1==0
                continue;
            end
            
            tmp2=PM2(tx-n,K2-i,K2-i,q);
            val=val+tmp1*tmp2;
        end
        
        sol=sol+val*binopdf(n,tx,g(1));
        
    end
    
    l2_prob(tx)=sol;
    
    disp(['Layer 2: ' num2str(tx) ' out of ' num2str(T_max)])
    
end

end

% Plotter for a nice graph!
function legend_h = plotter(l1_prob,l2_prob,T_min,T_step,T_max,K1,K2,g)

% Replace 0 with NaN in (l1_prob,l2_prob) for prettier plot
for k=1:length(l1_prob)
    
    if l1_prob(k)==0
        l1_prob(k)=NaN;
    end
end

for k=K2:length(l2_prob)
    if l2_prob(k)==0
        l2_prob(k)=NaN;
    end
    
end


% Plotting
figure(1)

if g(1,1)==0.3
    l1_style='-^';
    l2_style='--^'; 
end

if g(1,1)==0.4
    l1_style='-s';
    l2_style='--s'; 
end

if g(1,1)==0.5
    l1_style='-o';
    l2_style='--o'; 
end

if g(1,1)==0.6
    l1_style='-d';
    l2_style='--d'; 
end

if g(1,1)==0.7
    l1_style='-v';
    l2_style='--v'; 
end


% Legend fix
% plot(-1,-1,'-o',-1,-1,'--o',-1,-1,'-s',-1,-1,'--s',-1,-1,'--^',-1,-1,'--^');

plot(1:length(l1_prob),l1_prob,l1_style,'Color','k','MarkerSize',4)
% plot(1:1:length(l1_prob),l1_prob(1:1:end),l1_style,'Color','k','MarkerSize',4)

plot(1:length(l2_prob),l2_prob,l2_style,'Color','k','MarkerSize',4)
% plot(1:1:length(l2_prob),l2_prob(1:1:end),l2_style,'Color','k','MarkerSize',4)

% Plot annotation
grid('on')
pbaspect([2.5 1 1])
legend_h = legend('location','NorthWest');
xlabel('Total number of received packets [-]')
ylabel('Decoding probability [-]')
set(gca,'XTick',0:10:T_max)
set(gca,'YTick',0:0.1:1)
xlim([T_min T_max])
ylim([0 1])

end

% The New helper function
function P = PM2(m,n,r,q)

global PP_q1 PP_q2

s1=zeros(n+1,1);
s1(1)=1;

if q==2^1
    val=(PP_q1{n}^m)*s1;
    P=val(r+1);
end

if q==2^8
    val=(PP_q2{n}^m)*s1;    
    P=val(r+1);
end


end

% Returns a transition matrix of dim 'n+1' by 'n+1' with q param
function M = TransitionMatrix(n,q)
P1=zeros(n+1,n+1);
for i=1:length(P1)
    P1(i,i)=1/(q^(n-(i-1)));
    if i<n+1
        P1(i+1,i)=1-P1(i,i);
    end 
end

M=P1;

end















